Hands-on Exercise 10

Author

Dew Stella Chan

Overview

pacman::p_load(tmap, sf, DT, stplanr, tidyverse,sp, performance, reshape2, ggpubr, tidyverse)

Preparing the Flow Data

odbus <- read_csv("data/aspatial/origin_destination_bus_202210.csv")
glimpse(odbus)
Rows: 5,122,925
Columns: 7
$ YEAR_MONTH          <chr> "2022-10", "2022-10", "2022-10", "2022-10", "2022-…
$ DAY_TYPE            <chr> "WEEKDAY", "WEEKENDS/HOLIDAY", "WEEKENDS/HOLIDAY",…
$ TIME_PER_HOUR       <dbl> 10, 10, 7, 11, 16, 16, 20, 7, 7, 11, 11, 8, 11, 11…
$ PT_TYPE             <chr> "BUS", "BUS", "BUS", "BUS", "BUS", "BUS", "BUS", "…
$ ORIGIN_PT_CODE      <dbl> 65239, 65239, 23519, 52509, 54349, 54349, 43371, 8…
$ DESTINATION_PT_CODE <dbl> 65159, 65159, 23311, 42041, 53241, 53241, 14139, 9…
$ TOTAL_TRIPS         <dbl> 2, 1, 2, 1, 1, 4, 1, 3, 1, 5, 2, 5, 15, 40, 1, 1, …
odbus$ORIGIN_PT_CODE <- as.factor(odbus$ORIGIN_PT_CODE)
odbus$DESTINATION_PT_CODE <- as.factor(odbus$DESTINATION_PT_CODE) 
odbus6_9 <- odbus %>%
  filter(DAY_TYPE == "WEEKDAY") %>%
  filter(TIME_PER_HOUR >= 6 &
           TIME_PER_HOUR <= 9) %>%
  group_by(ORIGIN_PT_CODE,
           DESTINATION_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))
datatable(odbus6_9)
write_rds(odbus6_9, "data/rds/odbus6_9.rds")
odbus6_9 <- read_rds("data/rds/odbus6_9.rds")

Working with Geospatial Data

busstop <- st_read(dsn = "data/geospatial",
                   layer = "BusStop") %>%
  st_transform(crs = 3414)
Reading layer `BusStop' from data source 
  `C:\dewschan\ISSS622-GAA\Hands-on_Ex\Hands-On_Ex10\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 5159 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 3970.122 ymin: 26482.1 xmax: 48284.56 ymax: 52983.82
Projected CRS: SVY21
mpsz <- st_read(dsn = "data/geospatial",
                   layer = "MPSZ-2019") %>%
  st_transform(crs = 3414)
Reading layer `MPSZ-2019' from data source 
  `C:\dewschan\ISSS622-GAA\Hands-on_Ex\Hands-On_Ex10\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS:  WGS 84
mpsz
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21 / Singapore TM
First 10 features:
                 SUBZONE_N SUBZONE_C       PLN_AREA_N PLN_AREA_C       REGION_N
1              MARINA EAST    MESZ01      MARINA EAST         ME CENTRAL REGION
2         INSTITUTION HILL    RVSZ05     RIVER VALLEY         RV CENTRAL REGION
3           ROBERTSON QUAY    SRSZ01  SINGAPORE RIVER         SR CENTRAL REGION
4  JURONG ISLAND AND BUKOM    WISZ01  WESTERN ISLANDS         WI    WEST REGION
5             FORT CANNING    MUSZ02           MUSEUM         MU CENTRAL REGION
6         MARINA EAST (MP)    MPSZ05    MARINE PARADE         MP CENTRAL REGION
7                   SUDONG    WISZ03  WESTERN ISLANDS         WI    WEST REGION
8                  SEMAKAU    WISZ02  WESTERN ISLANDS         WI    WEST REGION
9           SOUTHERN GROUP    SISZ02 SOUTHERN ISLANDS         SI CENTRAL REGION
10                 SENTOSA    SISZ01 SOUTHERN ISLANDS         SI CENTRAL REGION
   REGION_C                       geometry
1        CR MULTIPOLYGON (((33222.98 29...
2        CR MULTIPOLYGON (((28481.45 30...
3        CR MULTIPOLYGON (((28087.34 30...
4        WR MULTIPOLYGON (((14557.7 304...
5        CR MULTIPOLYGON (((29542.53 31...
6        CR MULTIPOLYGON (((35279.55 30...
7        WR MULTIPOLYGON (((15772.59 21...
8        WR MULTIPOLYGON (((19843.41 21...
9        CR MULTIPOLYGON (((30870.53 22...
10       CR MULTIPOLYGON (((26879.04 26...
mpsz <- write_rds(mpsz, "data/rds/mpsz.rds")
busstop_mpsz <- st_intersection(busstop, mpsz) %>%
  select(BUS_STOP_N, SUBZONE_C) %>%
  st_drop_geometry()
datatable(busstop_mpsz)
od_data <- left_join(odbus6_9 , busstop_mpsz,
            by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>%
  rename(ORIGIN_BS = ORIGIN_PT_CODE,
         ORIGIN_SZ = SUBZONE_C,
         DESTIN_BS = DESTINATION_PT_CODE)
duplicate <- od_data %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()
od_data <- unique(od_data)
od_data <- left_join(od_data , busstop_mpsz,
            by = c("DESTIN_BS" = "BUS_STOP_N")) 
duplicate <- od_data %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()
od_data <- unique(od_data)
od_data <- od_data %>%
  rename(DESTIN_SZ = SUBZONE_C) %>%
  drop_na() %>%
  group_by(ORIGIN_SZ, DESTIN_SZ) %>%
  summarise(MORNING_PEAK = sum(TRIPS))
write_rds(od_data, "data/rds/od_data_fii.rds")
od_data_fii <- read_rds("data/rds/od_data.rds")

Visualising Spatial Interaction

od_data_fij <- od_data[od_data$ORIGIN_SZ!=od_data$DESTIN_SZ,]
write_rds(od_data_fij, "data/rds/od_data_fij.rds")
od_data_fij <- read_rds("data/rds/od_data_fij.rds")

Creating desire lines

flowLine <- od2line(flow = od_data_fij, 
                    zones = mpsz,
                    zone_code = "SUBZONE_C")
write_rds(flowLine, "data/rds/flowLine.rds")
flowLine <- read_rds("data/rds/flowLine.rds")
tm_shape(mpsz) +
  tm_polygons() +
flowLine %>%  
tm_shape() +
  tm_lines(lwd = "MORNING_PEAK",
           style = "quantile",
           scale = c(0.1, 1, 3, 5, 7, 10),
           n = 6,
           alpha = 0.3)

tm_shape(mpsz) +
  tm_polygons() +
flowLine %>%  
  filter(MORNING_PEAK >= 5000) %>%
tm_shape() +
  tm_lines(lwd = "MORNING_PEAK",
           style = "quantile",
           scale = c(0.1, 1, 3, 5, 7, 10),
           n = 6,
           alpha = 0.3)

#Hands Exercise 10b

Computing Distance Matrix

mpsz <- read_rds("data/rds/mpsz.rds")
mpsz
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21 / Singapore TM
First 10 features:
                 SUBZONE_N SUBZONE_C       PLN_AREA_N PLN_AREA_C       REGION_N
1              MARINA EAST    MESZ01      MARINA EAST         ME CENTRAL REGION
2         INSTITUTION HILL    RVSZ05     RIVER VALLEY         RV CENTRAL REGION
3           ROBERTSON QUAY    SRSZ01  SINGAPORE RIVER         SR CENTRAL REGION
4  JURONG ISLAND AND BUKOM    WISZ01  WESTERN ISLANDS         WI    WEST REGION
5             FORT CANNING    MUSZ02           MUSEUM         MU CENTRAL REGION
6         MARINA EAST (MP)    MPSZ05    MARINE PARADE         MP CENTRAL REGION
7                   SUDONG    WISZ03  WESTERN ISLANDS         WI    WEST REGION
8                  SEMAKAU    WISZ02  WESTERN ISLANDS         WI    WEST REGION
9           SOUTHERN GROUP    SISZ02 SOUTHERN ISLANDS         SI CENTRAL REGION
10                 SENTOSA    SISZ01 SOUTHERN ISLANDS         SI CENTRAL REGION
   REGION_C                       geometry
1        CR MULTIPOLYGON (((33222.98 29...
2        CR MULTIPOLYGON (((28481.45 30...
3        CR MULTIPOLYGON (((28087.34 30...
4        WR MULTIPOLYGON (((14557.7 304...
5        CR MULTIPOLYGON (((29542.53 31...
6        CR MULTIPOLYGON (((35279.55 30...
7        WR MULTIPOLYGON (((15772.59 21...
8        WR MULTIPOLYGON (((19843.41 21...
9        CR MULTIPOLYGON (((30870.53 22...
10       CR MULTIPOLYGON (((26879.04 26...

Converting from sf data.table to SpatialPolygonsDataFrame

mpsz_sp <- as(mpsz, "Spatial")
mpsz_sp
class       : SpatialPolygonsDataFrame 
features    : 332 
extent      : 2667.538, 56396.44, 15748.72, 50256.33  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 6
names       : SUBZONE_N, SUBZONE_C, PLN_AREA_N, PLN_AREA_C,       REGION_N, REGION_C 
min values  : ADMIRALTY,    AMSZ01, ANG MO KIO,         AM, CENTRAL REGION,       CR 
max values  :    YUNNAN,    YSSZ09,     YISHUN,         YS,    WEST REGION,       WR 
dist <- spDists(mpsz_sp, 
                longlat = FALSE)
head(dist, n=c(10, 10))
           [,1]       [,2]      [,3]      [,4]       [,5]      [,6]      [,7]
 [1,]     0.000  3926.0025  3939.108 20252.964  2989.9839  1431.330 19211.836
 [2,]  3926.003     0.0000   305.737 16513.865   951.8314  5254.066 16242.523
 [3,]  3939.108   305.7370     0.000 16412.062  1045.9088  5299.849 16026.146
 [4,] 20252.964 16513.8648 16412.062     0.000 17450.3044 21665.795  7229.017
 [5,]  2989.984   951.8314  1045.909 17450.304     0.0000  4303.232 17020.916
 [6,]  1431.330  5254.0664  5299.849 21665.795  4303.2323     0.000 20617.082
 [7,] 19211.836 16242.5230 16026.146  7229.017 17020.9161 20617.082     0.000
 [8,] 14960.942 12749.4101 12477.871 11284.279 13336.0421 16281.453  5606.082
 [9,]  7515.256  7934.8082  7649.776 18427.503  7801.6163  8403.896 14810.930
[10,]  6391.342  4975.0021  4669.295 15469.566  5226.8731  7707.091 13111.391
           [,8]      [,9]     [,10]
 [1,] 14960.942  7515.256  6391.342
 [2,] 12749.410  7934.808  4975.002
 [3,] 12477.871  7649.776  4669.295
 [4,] 11284.279 18427.503 15469.566
 [5,] 13336.042  7801.616  5226.873
 [6,] 16281.453  8403.896  7707.091
 [7,]  5606.082 14810.930 13111.391
 [8,]     0.000  9472.024  8575.490
 [9,]  9472.024     0.000  3780.800
[10,]  8575.490  3780.800     0.000
sz_names <- mpsz$SUBZONE_C
colnames(dist) <- paste0(sz_names)
rownames(dist) <- paste0(sz_names)
distPair <- melt(dist) %>%
  rename(dist = value)
head(distPair, 10)
     Var1   Var2      dist
1  MESZ01 MESZ01     0.000
2  RVSZ05 MESZ01  3926.003
3  SRSZ01 MESZ01  3939.108
4  WISZ01 MESZ01 20252.964
5  MUSZ02 MESZ01  2989.984
6  MPSZ05 MESZ01  1431.330
7  WISZ03 MESZ01 19211.836
8  WISZ02 MESZ01 14960.942
9  SISZ02 MESZ01  7515.256
10 SISZ01 MESZ01  6391.342

16.5.5 Updating intra-zonal distances

distPair %>%
  filter(dist > 0) %>%
  summary()
      Var1             Var2             dist        
 MESZ01 :   331   MESZ01 :   331   Min.   :  173.8  
 RVSZ05 :   331   RVSZ05 :   331   1st Qu.: 7149.5  
 SRSZ01 :   331   SRSZ01 :   331   Median :11890.0  
 WISZ01 :   331   WISZ01 :   331   Mean   :12229.4  
 MUSZ02 :   331   MUSZ02 :   331   3rd Qu.:16401.7  
 MPSZ05 :   331   MPSZ05 :   331   Max.   :49894.4  
 (Other):107906   (Other):107906                    
distPair$dist <- ifelse(distPair$dist == 0,
                        50, distPair$dist)
distPair %>%
  summary()
      Var1             Var2             dist      
 MESZ01 :   332   MESZ01 :   332   Min.   :   50  
 RVSZ05 :   332   RVSZ05 :   332   1st Qu.: 7097  
 SRSZ01 :   332   SRSZ01 :   332   Median :11864  
 WISZ01 :   332   WISZ01 :   332   Mean   :12193  
 MUSZ02 :   332   MUSZ02 :   332   3rd Qu.:16388  
 MPSZ05 :   332   MPSZ05 :   332   Max.   :49894  
 (Other):108232   (Other):108232                  
distPair <- distPair %>%
  rename(orig = Var1,
         dest = Var2)
write_rds(distPair, "data/rds/distPair.rds") 
distPair <- read_rds("data/rds/distPair.rds")

Preparing flow data

od_data_fii <- read_rds("data/rds/od_data_fii.rds")
flow_data <- od_data_fii %>%
  group_by(ORIGIN_SZ, DESTIN_SZ) %>% 
  summarize(TRIPS = sum(MORNING_PEAK))
head(flow_data, 10)
# A tibble: 10 × 3
# Groups:   ORIGIN_SZ [1]
   ORIGIN_SZ DESTIN_SZ TRIPS
   <chr>     <chr>     <dbl>
 1 AMSZ01    AMSZ01     1998
 2 AMSZ01    AMSZ02     8289
 3 AMSZ01    AMSZ03     8971
 4 AMSZ01    AMSZ04     2252
 5 AMSZ01    AMSZ05     6136
 6 AMSZ01    AMSZ06     2148
 7 AMSZ01    AMSZ07     1620
 8 AMSZ01    AMSZ08     1925
 9 AMSZ01    AMSZ09     1773
10 AMSZ01    AMSZ10       63
flow_data$FlowNoIntra <- ifelse(
  flow_data$ORIGIN_SZ == flow_data$DESTIN_SZ, 
  0, flow_data$TRIPS)
flow_data$offset <- ifelse(
  flow_data$ORIGIN_SZ == flow_data$DESTIN_SZ, 
  0.000001, 1)
flow_data$ORIGIN_SZ <- as.factor(flow_data$ORIGIN_SZ)
flow_data$DESTIN_SZ <- as.factor(flow_data$DESTIN_SZ)
flow_data1 <- flow_data %>%
  left_join (distPair,
             by = c("ORIGIN_SZ" = "orig",
                    "DESTIN_SZ" = "dest"))

Preparing Origin and Destination Attributes

pop <- read_csv("data/aspatial/pop.csv")

Geospatial data wrangling

pop <- pop %>%
  left_join(mpsz,
            by = c("PA" = "PLN_AREA_N",
                   "SZ" = "SUBZONE_N")) %>%
  select(1:6) %>%
  rename(SZ_NAME = SZ,
         SZ = SUBZONE_C)
flow_data1 <- flow_data1 %>%
  left_join(pop,
            by = c(ORIGIN_SZ = "SZ")) %>%
  rename(ORIGIN_AGE7_12 = AGE7_12,
         ORIGIN_AGE13_24 = AGE13_24,
         ORIGIN_AGE25_64 = AGE25_64) %>%
  select(-c(PA, SZ_NAME))
flow_data1 <- flow_data1 %>%
  left_join(pop,
            by = c(DESTIN_SZ = "SZ")) %>%
  rename(DESTIN_AGE7_12 = AGE7_12,
         DESTIN_AGE13_24 = AGE13_24,
         DESTIN_AGE25_64 = AGE25_64) %>%
  select(-c(PA, SZ_NAME))

Calibrating Spatial Interaction Models